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Abstract: 

High-resolution signal parameter estimation is a 
problem of significance in many signal processing 
applications. Such applications include direction- 
of-arrival (DOA) estimation, system identification, 
and time series analysis. A novel approach to the 
genera! Problem of signal parameter estimation is 
described. Although discussed in the context of 
direction-of-arrival estimation, FSPRIT can be 
applied to a wide variety of problems including 
accurate detection and estimation of sinusoids in 
noise. It exploits an underlying rotational 
invariance among signal subspaces induced by an 
array of sensors with a rotational invariance 
structure. The technique, when applicable, mani- 
fests significant performance and computational 
advantages over previous algorithms such as 
MEM, Capon's MLM, and MI SIC. 
Keywords: DOA, FSPRIT, MEM, MLM 

1. INTRODUCTION 

In many practical signal processing problems, 
the objective is to estímate from measurements a set 
of constant parameters upon which the received 
signáis depend. For example, high-resolution 
direction-of-arrival (DOA) estimation is important in 
many sensor systems such as radar, sonar, electronic 
surveillance, and seismic exploration. High-resolution 
frequency estimation is important in numerous 
applications, recent examples of which include the 
design and control of robots and large flexible space 
structures. In such problems, the functional form of 
the underlying signáis can often be assumed to be 
known (e.g. narrow-band plañe waves, cisoids). The 
quantities to be estimated are parameters (e.g., 
frequencies and DOA's of plañe waves, cisoid 
frequencies) upon which the sensor outputs depend, 
and these parameters are assumed to be constant. 1 

There have been several approaches to such 
problems including the so-called máximum likelihood 
(ML) method of Capón (1969) and Burg's (1967) 
máximum entropy (ME) method. Although often 
successful and widely used, these methods have 
certain fundamental limitations (especially bias and 
sensitivity in parameter estimates), largely because 
they use an incorrect model (eg. AR rather than 
special ARMA) of the measurements. Pisa renko 
(1973) was one of the first to exploit the structure of 



the data model, doing so in the context of estimation 
of parameters of cisoids in additive noise using a 
covariance approach. Schmidt (1977) and 
independently Bienvenu (1979) were the first to 

correctly exploit the measurement model in the case 
of sensor arrays of arbitrary form. Schmidt, in 
particular, accomplished this by first deriving a 
complete geometric solution in the absence of noise, 
then cleverly extending the geometric concepts to 
obtain a reasonable approximate solution in the 
presence of noise. The resulting algorithm was called 
MUSIC (Múltiple Signal Classification) and has been 
widely studied. In a detailed evaluation based on 
thousands of simulations, M.I.T/s Lincoln Laboratory 
concluded that, among currently accepted high- 
resolution algorithms, MUSIC was the most 
promising and a leading candidate for further study 
and actual hardware implementation. However, 
although the performance advantages of MUSIC are 
substantial, they are achieved at a considerable cost 
in computation (searching over parameter space) and 
storage (of array calibration data). 

In this paper, a new algorithm (ESPRIT) that 
dramatically reduces these computation and storage 
costs is presented. In the context of DOA estimation, 
the reductions are achieved by requiring that the 
sensor array possess a displacement invariance, i.e., 
sensors occur in matched pairs with identical 
displacement vectors. Fortunately, there are many 
practical problems in which these conditions are or 
can be satisfied. In addition to obtaining signal 
parameter estimates efficiently, optimal signal copy 
vectors for reconstructing the signáis are elements of 
the ESPRIT solution as well. ESPRIT is also 
manifestly more robust (i.e., less sensitive) with 
respect to array imperfections than previous 
techniques including MUSIC. 

To make the presentation as clear as possible, 
an attempt is made to adhere to a somewhat standard 
notational convention. Lowercase boldface italic 
characters will generally refer to vectors. Uppercase 
boldface italic characters will generally refer to 
matrices. For either real- or complex-valued 
matrices., (•)* will be used to denote the Hermitian 
conjúgate (or complex-conjugate transpose) 
operation. Eigenvalues of square Hermitian matrices 
are assumed to be ordered in decreasing magnitude, 
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as are "the singular valúes of nonsquare matrices. 
Knowledge of the fundamental theorems of matrix 
algebra dealing with ei-gendecompositions and 
singular valué decompositions (SVD) is assumed (cf. 
[2]). 

THE DATA MODEL 

Although FSPRIT is generally applicable to a 
wide variety of problems, for illustrative purposes the 
discussions herein focus on DOA estimation. In many 
practical signal processing applications, data from an 
array of sensors are collected, and the objective is to 
lócate point sources assumed to be radiating energy 
that is detectable by the sensors (cf. Fig. 1). 
Mathematically, such problems are quite simply, 
although abstractly, modeled using Green's functions 
for the particular differential operator that describes 
the physics of radiation propagation from the sources 
to the sensors. For the intended applications, 
however, a few reasonable assumptions can be 
invoked to make the problem analytically tractable. 

The transmission médium is assumed to be 
isotropic and nondispersive so that the radiation 
propagates in straight lines, and the sources are 
assumed to be in the far-field of the array. 
Consequently, the radiation impinging on the array is 
in the form of a sum of plañe waves. For simplicity, it 
will initially be assumed that the problem is planar, 
thus reducing the location parameter space to a 
single-dimensional subset of (R, Le., 9¡e [- n, n], 
where 9¡ is the direction-of-arrival (DOA) of the ith 
source. 

The signáis are assumed to be narrow-band 
processes, and can be considered to be sample 
functions of a stationary stochastic process or 
deterministic functions of time. This can represent a 
significant difference when estimating the signáis 
themselves, leading to entirely different estimation 
algorithms. However, as far as estimation of signal 
parameters such as DOA is concerned, both as- 
sumptions lead to the same (suboptimal) algorithm 
under certain conditions (e.g., persistent excitation). 

Since the narrow-band signáis are assumed to 
have the same known center frequency 2 (ff> ), the ith 
signal can be written as s¡ (t) = u¡ (t) eos (a>ot + v¡(t), 
where u¡ (t) and v¡ (t) are slowly varying 3 functions 
of time that define the amplitude and phase of s¡ (t), 
respectively. For such signáis, it is often more 
convenient to use the complex envelope 
representation [3] in which s,(t) = Re (s(t)}, where 
s(t) = u(t) e I(to0t + vi(t) . Noting that the narrowband 
assumption implies u(t) = u(t - x) and v(t)= v(t - x) for 
all possible propagation delays x, the effect of a time 
delay on the received waveforms is simply a phase 
shift, i.e., s(t - x)= s(t) e" Jt0OT . The result is that x k (t), 
the complex signal output of the kth sensor at time t, 
can be written as 



X k (t) = ¿a k (6 1 )s i (t-i 

¡=i 

= ¿a k( e l )s l (t )e - J - (e ' ) 
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where x k (9) is the propagation delay between a 
reference point and the kth sensor for the ith 
wavefront impinging on the array from direction 
9,a k (9i) is the corresponding sensor element complex 
response (gain and phase) at frequency (0 o , and there 
are assumed to be d point sources present. 

Employing vector notation for the outputs of the m 
sensors, the data model becomes 



d 



(2) 



jco0Tm(9i) -,T 



y (3) 



x(t) = 2,a(9 i )s i (t) 

¡=i 

Where 

a(0 1 ) = [a 1 (0 1 )e- J ^ 1(ei) ),...,a m (9,)e- 

often termed the array response or array steering 
vector for direction 9j. Setting A(9) = [a(9 1)j 
a (6d)L s(t) = [s i(t), s d (t)] T , and adding measurement 
noise n(t), the measurement model for the passive 
sensor array narrow-band signal processing problem 
is 



x(t) =A(9) s(t) + n(t). 



Note that x(t), n(t)e C m , s(t) e c d 
and it will be assumed that m > d. 



III. 



(4) 



and A(9)e C n 



THE GEOMETRIC APPROACH 

In 1977, Schmidt [4] developed the MUSIC 
(Múltiple Signal Classification) algorithm by taking a 
geometric view of the signal parameter estimation 
problem. One of the major breakthroughs afforded by 
the MUSIC algorithm was the ability to handle 
arbitrary arrays of sensors. 




Until the mid-1970's. direction finding 
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tecTTiiiques required knowledge of the array 
directional sensitivity pattern in analytical form, and 
the task of the antenna designer was to build an array 
of antennas with a pre-specified sensitivity pattern. 
The work of Schmidt essentially relieved the designer 
from such constraints by exploiting the reduction in 
analytical complexity that could be achieved by 
calibrating the array. Thus, the highly nonlinear 
problem of calculating the array response to a signal 
from a given direction was reduced to that of measur- 
ing and storing the response. Although MUSIC did 
not mitigate the computational complexity of solution 
to the DOA estimation problem, it did extend the 
applicability of high-resolution DOA estimation to 
arbitrary arrays of sensors. 

A. Array Manifolds and Signal Subspaces 

To introduce the concepts of the array 
manifold and the signal subspace, recall the noise 
free data model x(t) = A(8) s(t). The vectors a (9,) 
eC". the columns of A (9) are elements of a set (not 
a subspace), termed the array manifold (a), composed 
of all array response (steering) vectors obtained as 9 
ranges over the entire parameter space. a is 
completely determined by the sensor directivity 
patterns and the array geometry, and can sometimes 
be computed analytically. However, for complex 
arrays that defy analytical description, a can be 
obtained by calibration (i.e., physical measurements). 

For azimuth-only DOA estimation, the array 
manifold is a one-parameter manifold that can be 
viewed as a rope weaving through C". For azimuth 
and elevation DOA estimation, the manifold is a sheet 
in C" To avoid ambiguities, it is necessary to assume 
that the map from 9 

¡il , 9d} to 3í {A(9)}, the subspace spanned by 
the columns of A (9), is one-to-one. This property can 
be ensured by proper array design. 



The key observation is that if x(t) = a (9) Se 
(t) is an appropriate data model (in the absence of 
noise) for a single signal, the data are confined to a 
one-dimensional sub-space of characterized by the 
vector a(9). For d signáis, the observed data vectors 
x(t) = A(9) s(t) are constrained to the d-dimensional 
subspace of C", termed the signal subspace (S x ), that 
is spanned by the d vectors a(9¡), the columns of 
A(9). 

B. Intersections as Solutions 

The concepts of an observed signal subspace 
and a calibrated array manifold permit an immediate 
visualization of the solution. In the absence of noise, 
the outputs of the sensor array lie in a d-dimensional 
subspace of C" the signal subspace (S x ) spanned by 
the columns of A(9). Once d independent vectors 
have been observed, S x is known, and intersections 
between the observed subspace and the array 
manifold yield the set of vectors from the array 
manifold that span the observed signal subspace. A 
three-sensor two-source example is graphically 
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depicted in Fig. 2. Assuming that the sensor array has 
been designed such that the map from parameters to 
array manifold vectors is unique, the parameters are 
immediately determined. 

Problems arise when only noisy 
measurements x(t) = A(9) s(t) + n(t) of the array 
output are available, since S x must be estimated- 
Imposing the constraint that the estimate S x be 
spanned by elements from a and assuming unknown 
deterministic signáis and Gaussian noise, a 
maximum-likelihood (ML) estimator can be 
formulated as described in Appendix A. However, the 
ML solution (of obtaining a set of vectors from the 
array manifold that best fits all the measurements) is 
computationally prohibitive in most practical 
applications. 

Schmidt's idea was to employ a suboptimal 
two-step procedure instead. First, an unconstrained 
set of d vectors that best fits all the measurements is 
found. Then points of closest approach of the space 
spanned by those vectors to the array manifold are 
sought. This procedure, although clearly suboptimal. 
retains some of the key properties of the ML solution, 
including the fact that the exact answer is obtained 
asymptotically as the number of measurements goes 
to infinity.) 

C. Estimating the Signal Subspace 

To obtain an unconstrained estimate of the 
signal sub-space, the least-squares (LS) criterion is 
most often employed. The idea is to find a set of d 
vectors that span a subspace of C" that best fits, in an 
LS sense, the observed data. Assuming the signáis 
and noise are zero-mean, a method can be derived by 
first examining the covariance matrix of the 
measurements. If the signáis are modeled as 
stationary (zero-mean) stochastic processes, they are 
assumed to be uncorrelated with the noise and possess 
a positive definite covariance matrix R ss > 0. If, on 
the other hand, a deterministic (zero-mean) signal 
model is chosen, a persistent excitation condition is 
imposed, i.e., 

N 



R, 



def | _ 

= lim gs(t) 



N 



N 



s*(t) 



is assumed to exist and be positive definite. 

Under the conditions given above, the 
covariance matrix of the measurements is given by 

R xx = E{xx*} = AR SS A* + a 2 Z n (5) 

The objective is to find a set of d linearly 
independent vectors that is contained in S x = R{A}, 
the subspace of the eigenvectors have been proposed 
by several investigators, e.g., Johnson 16]. 

Although conceptually simple, the one- 
dimensional MUSIC measure has several drawbacks. 
Primarily, problems in the finite measurement case 
arise from the fact that since d signáis are known to 
be present, d parameter estimates. {9|...,9 d } 
should be sought simultaneously by maximizing an 
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apprópriate functional rather than obtaining estimates 
one at a time as is done in the search over Pm(8)- 
However, multidimensional searches are 
exponentially more expensive than one-dimensional 
searches. The price paid for the computational 
reduction achieved by employing a one-dimensional 
search for d parameters is that the method is finite- 
sample-biased in the múltiple source environment (cf. 
[1]). Furthermore, in low SNR scenarios and in 
situations where even small sensor array errors are 
present, the ability of the conventional MUSIC 
spectrum to resolve closely spaced sources (i.e., 
observe múltiple peaks in the measure) is severely 
degraded (cf. [1]). Nevertheless, it should be 
emphasized that in spite of these drawbacks, MUSIC 
has been shown to outperform previous techniques 
(cf. [7]). 

Finally, as indicated above, MUSIC 
asymptotically yields unbiased parameter estimates 
since as the amount of data becomes infinite, errors in 
the estimate S x vanish If the noise is spatially 
Gaussian and temporally independent, the distribution 
of the eigenvectors of the sample covariance is 
asymptotically Gaussian, with mean equal to the true 
eigenvectors (assuming distinct eigenvalues) and 
covariance that go to zero [9], [10]. Thus, the 
estimated signal subspace converges in mean-square 
to the true signal subspace, and the parameter es- 
timates converge to the true valúes as well. 

E. Summary of the MUSIC Algorithm 

The following is a summary of the MUSIC 
algorithm based on the covariance approach described 
above. 



1) Collect the data and estimate 



R* 



def 

_ E{xx*j 



AR SS A* + ct % 



2) 



3) 
4) 



denoting the estimate R xx . 

Solve for the eigensystem; RxxE = £„E /l , where 
is A = diag • • • X m }, %i ... X m , and E = [ei.. 
I... lej. 

Estimate the number of sources d. 
Evalúate 



Pm(0)= 



a* O) a (0) 



a*(0)E N E N a(0) 



where E N = £ n [ e d+1 l ... - le m ]. 

5) Find the d (largest) peaks of P(9) to obtain estimates of 
the parameters. 

An alternative to first forming the 
measurement covariance matrix and then performing 
an eigendecom posit ion is to opérate directly on the 
measurements using the singular valué decomposition 
(SVD). In addition to avoiding squaring the data, this 
approach has a nice geometric interpretation. Letting 
X eC mxN denote the data matrix, the objective is to 
obtain a set of vectors spanning the column space of 
the rank d matrix, X eC mxN , that best approximates X 
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in a least-squares sense. The solution is given by the 
d left singular vectors of X corresponding to the d 
largest singular valúes. 12 

It is easily demonstrated that the 
eigendecomposition and the SVD yield the same 

subspace estimate. If the SVD of X/VÑ" is given by 
UEV*, 

1 



N 



XX* = U£ Z U* = Rxx. 



since £ is diagonal and real, and U and V are unitary. 
Thus, the left singular vectors, U, of X are the right 
eigenvectors of R xx , the sample covariance matrix. 
Thus, in the absence of finite precisión effects in 
computing the decompositions, the subspace 
estimates from both techniques are identical. 

There are, however, significant computational 
differences. The computation of the full SVD of X is 
of order mN 2 which can be significantly larger than 
the 0(m 3 ) operations required for an 
eigendecomposition of Rxx. The increase in 
computation is due to the fact that the full SVD is 
also obtaining a set of d vectors, the first d columns 
of V 6 N * N that span the d-imensional subspace of N- 
dimensional space spanned by the d signal vectors, 
vectors in C N whose components are the samples of 
the underlying signáis. 13 If the information of the full 
SVD is not required, partial SVD algorithms that 
compute only the left singular vectors and singular 
valúes can be cm-ployed resulting in substantial 
computational savings [11]. 

If N is the number of measurements to be 
processed, and m is the number of elements in each 
sample vector, forming the sample covariance matrix 
requires on the order of Nm 2 operations. 
Eigendecompositions of matrices in 5í mxm require on 
the order of 10m 3 operations, whereas a standard 
SVD of X e9í mxN requires * 2Nm 2 + 4w 3 (cf. [2, p. 
175]) if only the singular valúes and left singular 
vectors are required. The computational effort is 
therefore on the same order for both the SVD and 
eigenvector approaches to reducing the measurements 
to the statistic E s used by ESPRIT and MUSIC Again, 
the major advantage of the SVD over the 
eigendecomposition is that the measurements are 
processed directly without squaring them. Thus, 
numerical problems associated with ill-conditioned 
matrices are mitigated to some extent by using the 
SVD. 

IV. ESPRIT 

Although MUSIC was the first of the high- 
resolution algorithms to correctly exploit the 
underlying data model of narrow-band signáis in 
additive noise, the algorithm has several limitations 
including the fact that complete knowledge of the 
array manifold is required, and that the search over 
parameter space is computationally very expensive. 
In this section, an approach (ESPRIT) to the signal 
parameter estimation problem that exploits sensor ar- 
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fay~in variances is described. ESPRIT is similar to 
MUSIC in that it correctly exploits the underlying 
data model, while manifesting significant advantages 
over MUSIC as described in Section I. 

To simplify the description of the basic ideas 
behind ESPRIT, much of the ensuing discussion is 
couched in terms of the problem of múltiple source 
direction-of-arrival (DOA) estimation from data 
collected by an array of sensors. For simplicity, 
discussions deal only with single-dimensional 
parameter spaces, e.g., azimuth-only direction finding 
(DF) of far-field point sources, since the basic 
concepts are most easily visualized in such spaces. 
Narrow-band signáis of known center frequency will 
be assumed. Recall that a DOA/DF problem is 
classified as narrow-band if the signal bandwidth is 
small compared to the inverse of the transit time of a 
wave front across the array, and the array response is 
not a function of frequency over the signal 
bandwidth. The generality of the fundamental 
concepts on which ESPRIT is based makes the 
extensión to higher spatial dimensions and to signáis 
containing múltiple frequencies possible. 

A. Array Geometry 

ESPRIT retains most of the essential features 
of the arbitrary array of sensors, but achieves a 
significant reduction in computational complexity by 
imposing a constraint on the structure of the sensor 
array, a constraint most easily described by an 
example. Consider a planar array of arbitrary 
geometry composed of m sensor doublets as shown in 
Fig. 3. The elements in each doublet have identical 
sensitivity patterns and are translationally separated 
by a known constant displacement vector A. Other 
than the obvious requirement that each sensor have 
nonzero sensitivity in all directions of interest, the 
gain, phase, and polarization sensitivity of the 
elements in the doublet are arbitrary. Furthermore, 
there is no requirement that any of the doublets 
possess the same sensitivity patterns although, as 
discussed in [1] and [14], there are advantages to 
employing arrays with such characteristics. 

B. The Data Model 

Assume that there are d< m narrow-band 
sources centered at frequency co and that the sources 
are located 
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Fig. 3. Sensor array geometry for múltiple source 
DOA estimation using ESPRIT. 

sufficiently far from the array such that in 
homogeneous isotropic transmission media, the 
wavefronts impinging on the array are planar, As 
before, the sources may be assumed to be stationary 
zero-mean random processes or deterministic signáis. 
Additive noise is present at all 2 m sensors and is 
assumed to be a stationary zero-mean random process 

with a spatial covariance a 

To describe mathematically the effect of the 
translational invariance of the sensor array, it is 
convenient to describe the array as being comprised 
of two subarrays, Z x and Z Y , identical in every 
respect although physically displaced (not rotated) 
from each other by a known displacement vector A of 
magnitude A. The signáis received at the ith doublet 
can then be expressed as 

X 1 (t) = ¿s k (t)a 1 (e k ) + n xl (t), 



d 



y,(t) = I> k (t)e 



jcü()Asin9k/c 



a,(9 k ) + n yi (t), 



k=l 



where e k is now the direction-of-arrival of the 
kth source relative to the direction of the translational 
displacement vector A. 

Since the sensor gain and phase patterns are arbitrary 
and since ESPRIT does not require any knowledge of 
the sensitivities, the subarray displacement vector A 
sets not only the scale for the problem, but the 
reference direction as well. The DOA estimates 
obtained are angles-of-arrival with respect to the 
direction of the vector A. A natural consequence of 
this fact is the necessity for a corresponding 
displacement vector for each dimensión in which 
parameter estimates are desired. 

Combining the outputs of each of the sensors 
in the two subarrays. the received data vectors can be 
written as follows: 

x(t)= As(t) + n x (t), (9) 

y(t) = A Os (t) + n, (r), (10) 

where the vector s(f) is the d x 1 vector of impinging 
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signáis (wavefronts) as observed at the reference 
sensor of subarray Z x . The signáis can be correlated 
in the sense E {s¡ (t) s ¡ (t)} # for i # j although the 
case of coherent (or fully correlated) sources is not 
considered herein (cf. [I, Section 7.11] for further 
discussion). The matrix d> is a diagonal d x d matrix 
of the phase delays between the doublet sensors for 
the d wavefronts, and is given by 



cD = diag {e 



jyi 



Jyd. 



(11) 



where y k = coA sin e k /c.cl> is a unitary matrix 
(operator) that relates the measurements from 
subarray Z x to those from subarray Z Y In the complex 
field, O is a simple scaling operator. However, it is 
isomorphic to the real two-dimensional rotation 
operator and is herein referred to as a rotation' 
operator. The unitary nature of O is a consequence of 
the narrow-band planewave assumption, an 
assumption that leads to unit-modulus eisoidal signáis 
in the spatial domain. In time series analysis, the 
diagonal elements of <B are potentially arbitrary 
complex numbers in which case <£ could be an 
expansive or contractive operator. 

Defining the total array output vector as Z(t), 
the sub-array outputs can be combined to yield 



z ( t ) = [ x y S]-Ás(t) + nz(t), 
Á = [U n 1 (t)=bíj] 



(12) 
(13) 



It is the structure of Á that is exploited to 
obtain estimates of the diagonal elements of <J> 
without having to know A. 

From (12), it is easily seen that the estimation 
problem posed is scale-invariant in the sense that 
absolute signal powers are not observable. For any 
nonsingular diagonal matrix, D, the data model is 
invariant with respect to the transformations s(t)— > D" 
1 s(t) and Á — > ÁD. Thus, estimates of the signáis and 
the associated array manifold vectors derived herein 
are to be interpreted modulo an arbitrary scale factor 
unless knowledge of the gain pattern of one of the 
sensors is available 

C. ESPRIT- The Invariance Approach 

The basic idea behind ESPRIT is exploit the 
rotational invariance of the underlying signal 
subspace induced by the translational invariance of 
the sensor array. The relevant signal subspace is the 
one that contains the outputs from the two subarrays 
described above, Z x and Z Y . Simultaneous sampling 
of the output of the arrays leads to two sets of 
vectors, E x and E Y , that span the same signal 
subspace (ideally, that spanned by the columns of A). 

The ESPRIT algorithm is based on the 
following results for the case in which the underlying 
2m-dimensional signal subspace containing the entire 
array output is known in the absence of noise, the 
signal subspace can be obtained as before by 
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collecting a sufficient number of measurements and 
finding any set of d linearly independent 
measurement vectors. These vectors span the d- 
dimensional subspace of C m spanned by Á. The signal 
subspace can also be obtained from knowledge of the 
covarianee of the measurements R zz = ÁR SS A* a 2 Z n . 
As discussed in Section III, if d < m (an assumption 
required in later developments), the 2 m - d smallest 
GE's of (Rzz,Z n ) are equal to a 2 . The GEV's 
corresponding to the d largest GE's are used to obtain 

E s = SJeJ.... Ie d l, where 5Í { E s } = 9Í { A }, the signal 
subspace. 

Since Jí {E s } = Jí {A }, there must exist a unique (re- 
call d<m), non singular T such that 



E s = AT 



(14) 



Furthermore, the invariance structure of the array 
implies E s can be decomposed into E x € C m * d and E Y 
£ £ M * D ( c f 2^ an( j z y subarrays) such that 







"AT 






A())T 



E s = „ A = . ,„ (15) 



from which it is easily seen that 

9t {E x }=5?{E y ) = 9í{A}. 

Since E x and E Y share a common column space, the 
rank of 

E XY det = [E x I E Y ] (16) 
is d, which implies there exists a unique 
(recall d < m ) rank d matrix 18 FGC 2dxd such that 

= [E X IE Y I F = E X F X + E Y Fy, (17) 
= A TF X + A<DTF Y (18) 
spans the null-space of [ E x \ E Y \. Defining 

def 

V=FxtF y ]-' (19) 
equation (18) can be rearranged to yield 

ATv|/= AÍ>T => ATv|/T 1 = A<D. 
Assuming A to be full rank implies 

Tv|/T"'=<D (21) 
Therefore, the eigenvalues of y must be equal 
to the diagonal elements of <D, and the columns of T 
are the eigenvectors of \|/ This is the key relationship 
in the development of ESPRIT and its properties. The 
signal parameters are obtained as nonlinear functions 
of the eigen-values of the operator \|/ that maps 
(rotates) one set of vectors [Ex) that span an m- 
dimensional signal subspace into another ( Ey). 

D. Estimating the Subspace Rotation Operator 

In practical situations where only a finite 
number of noisy measurements are available, E s 



estimated from the covarianee matrices of the 
measurements R zz or, equivalently, from the data 
matrix Z. The result is that 9í{E s ) is only an estimate 

of S z , and with probability one, 9Í {E s } # 5í { A }. 
Furthermore, 9í{E x ) # 9?{E Y }. Thus, the objective of 
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nndíñg a y such that E x y = E y is no longer 
achievable. A criterion for obtaining a suitable 
estímate must be formulated. The most commonly 
employed criterion for problems of this nature is the 
least-squares (LS) criterion. 

The standard LS criterion applied to the 
model AX = B to obtain an estimate of X assumes A 
is known and the error is to be attributed to B. 
Assuming the set of equations is overdetermined, the 
columns of A are linearly independent, and the noise 
in the elements of B is zero-mean and E{b¡j b k] } = 
CT 2 5 ik 5ji the LS solution is 

"X = IA*AI _1 A*B 

It is easily verified that the estimate is 
unbiased and minimum variance. The extensión to 
arbitrary, but known, covariance of the rows of B is 
straightforward and leads to the weighted least- 
squares (WLS) solution. If both A and B are noisy, 
however, the LS solution is known to be biased. 

Since it is not difficult to argüe that the 
estimates E x and E y are equally noisy, the LS criterion 
is clearly in appropriate. A criterion that takes into 
account noise on both A and B is the total least- 
squares (TLS) criterion. The TLS criterion can be 
stated [21 as finding residual matrices R A and R B of 
minimum Frobenius norm, and "X such that 

IA + R A I X = B + R B 

(22) 

This criterion is easily shown to be equivalent 
to replacing the zero matrix in (17) by a matrix of 
errors the Frobenius norm of which is to be 
minimized (i.e., total least-squared error). If the 
covariance of the errors, specifically the rows of [Ra 
IR B ] is known to withm a scale factor, the TLS 
estimate of X is strongly consistent [8]. 

Appending a nontriviality constraint F*F =1 
to eliminate the zero solution and applying standard 
Lagrange techniques leads to a solution for F given 
by the eigenvectors corresponding to the d smallest 
eigenvalues of E* X yE xy . The eigenvalues of y as 
defined above and cal culated from the estimates F x 
and F y are taken as estimates of the diagonal elements 
of (D 

E. Summary of the TLS ESPRIT Covariance 

Algorithm 

The TLS ESPRIT algorithm based on a 
covariance formulation can be summarized as 
follows. 

1) Obtain an estimate of R zz , denoted "R zz , from 
the measurements Z 

2) Compute the generalized eigendecomposition 
of{R zz , 2 n } 



Rzz E - 2 n E a 

where A = diag { Xi,... X 2m ], X\> ... > X 2m , and E = 
leil...le 2m l. 

3) Estimate the number of sources "d. 
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4) Obtain the signal subspace estimate "S z =3í{E s j 
and decompose it to obtain E x and E y , where 



def 



5) Compute the eigendecomposition (X { > . . . .> X 2á ) 
E* 



def 

ExY^XY = 



[E 



x IE Y ] = EAE : 



and Partition E into "d x"d submatrices. 



def 

E = 



^11 ^12 

E21 E 22 



6) Calcúlate the eigenvalues of v|/= -E 12 



v 22 



k ); > k = X k (-E 12 E 22 ) _1 for all k 
for DOA estimation, 9 k = sin" 1 {c arg 



7) Estimate "9 k : 
= 1,2 "d e.g. 
(4 k }/(co A)}. 

For arrays with múltiple in variances, such as 
uniform linear arrays, the decomposition of E s into Ex 
and E Y is not unique. See [141 and [17] for more 
details concerning múltiple invanance ESPRIT. 

In many instances, it is preferable to avoid 
forming covariance matrices, and instead to opérate 
directly on the data as discussed in Section 111. This 
approach leads to (generalized) singular valué 
decompositions (GSVD's) of data matrices, and a 
GSVD variant of ESPRIT discussed in detail in [ 1]. 

From the key relation (21). several other 
quite striking results can be derived. For example, not 
only is knowledge of the array manifold not required, 
but the elements thereof associated with the estimated 
signal parameters (DOA's) can be estimated if 
desired. The same is true of the source correlation 
matrix, knowledge of which is not needed in ESPRIT. 

E. Array Calibration 

Using the TLS formulation of ESPRIT, the 
array manifold vectors associated with each signal 
(parameter) can be estimated (to within an arbitrary 
scale factor). From (21), the right eigenvectors of v|/ 
are given by E V =T" . This result can be used to obtain 
estimates of the array manifold vectors as 

E S E ¥ = ATT -1 =A 
(23) 

No assumption concerning the source 
covariance is required. 

Although simple to compute, this estimate 
will not in general conform to the invariance structure 
of the array in the presence of noise. In low SNR 
scenarios, the deviation from the assumed structure 

A=[A T l(A<D) T ] T may be significant. In such 
situations, improved estimates of the array manifold 
vectors can be obtained by employing the formulation 
discussed in [18]. 

G. Signal Copy 
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~~ ln many practical applications, not only the 
signal parameters, but the signáis themselves, are of 
interest. Estimation of the signáis as a function of 
time from an estimated DOA is termed signal copy. 
The basic objective is to obtain estimates "s(t) of the 

signáis s(t) from the array output, z(t) = A s(t) + n(t). 
Employing a linear estimator, a squared-error cost 
criterion in the metric of the noise (which is ML if 
the noise is Gaussian), and conditioning on 

knowledge of A , leads to the estimate "s(t), the 
vector of coefficients resulting from the oblique pro- 
jection of z(t) onto the space spanned by the columns 
of A (cf. the Appendix). The resulting weight matrix 
W (i.e., the linear estimator) whose ith column is a 
weight vector that can be used to obtain an estimate 
of the signal from the itha estimated DOA and reject 
those from the other 

DOA's is given by 



w=v _1 a A*y~'Á 



(24) 

In terms of quantities already available, (24) can be 
written as 



W= y-' E [ E *y-' E VV-* 



(25) 

using (23) to estimate A. This equivalence is easily 
established since from (21) it follows that the right 
eigenvectors of \\i equal T . Combining this fact with 
E s =ÁT and substituting in (25) yields 



w*=e: 



e*y -1 e l" 1 E*y _1 



(26) 



Note that the optimal copy vector is a vector 
that is Z _1 n orthogonal to all but one of the vectors in 

the columns of A since W* A = I 

There is, of course, a total least-squares 

alternative to conditioning on knowledge of A . Since 

only estimates of A are available, in low SNR 
scenarios where accurate signal estimates are desired, 
the TLS approach yields improved estimates at the 
cost of increased computation. Although not derived 
herein, "s{t) can be obtained by performing a 
(generalized) singular valué decomposition of 

[A\z(t)l.] The right singular vector corresponding to 
the smallest singular valué yields ~s(t) as the first d 
elements after normalizing the last element to 
unity."" 

H. Source Correlation Estimation 

There are several approaches that can be used 
to estimate the source correlations. The most 
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straightforward is to simply note that the optimal 
signal copy matrix W obtained above removes the 
spatial correlation in the observed measurements [cf. 
(25)]. Thus, W*C ZZ W = DSD* where S is the source 
correlation (not covariance) matrix, C zz = R zz -cr 2 Z n , 
and the diagonal factor D accounts for arbitrary 
normalization of the columns of W. Note that when 
R zz must be estimated, a manifestly rank d estimate 



^ZZ 

A ( s d) 



can be used [cf. (6)], where 



diag {^i,... X¿, }_and X 17 is a GE of (R zz , Z n )- 
Combining this with E s = A T gives 



DSD* = T[ A ( s d) -CT2Id)T*. 



(27) 



If a gain pattern for one of the elements is 
known, specifically if the gain I a¡(9 k )l is known for 
all 9 K associated with sources whose power is to be 
estimated, then source power estimation is possible 
since the array manifold vectors can now be obtained 
with proper scaling. 

V. SlMULATlON Results 

Many simulations have been conducted 
exploring different aspects of ESPRIT and making 
comparisons to other techniques (cf. [1]). Herein, 
only one of the scenarios, but one that addresses 
several issues that arise in a practical implementation 
of ESPRIT, is presented. Thus, sensor gain and phase 
errors, as well as sensor spacing errors, are included. 
Furthermore, unequal source powers and a high 
degree of source correlation are assumed. 

More specifically, the array chosen was a ten- 
element array with doublet spacing X/4 and the five 
doublets randomly spaced on a line resulting in an 
aperture of approximately 4-X. Two sources were 
located at 24° and 28° (approximately 0.3 Rayleigh or 
3 dB beamwidth separation), and were of unequal 
powers, 20 dB and 15 dB, respectively. Sensor errors 
were introduced by zero-mean normal random 
additive errors with sigmas of 0.1 dB in amplitude 
and 2 o in phase 23 (independent of angle). Sensor 
location errors (along the axis of the array) with 
sigma 0.005 (X/2) were included as well. The sources 
were 90 percent temporally correlated and 5000 triáis 
were run. A histogram of the results is given in Fig. 
4. 

The number of sources was assumed to be 
known in the implementation of both MUSIC and 
ESPRIT. The indicated failure rate for MUSIC of 37 
percent is the percentage of triáis in which the 
conventional MUSIC spectrum did not exhibit two 
peaks in the interval [20°, 32° J. This, of course, is 
not an issue in ESPRIT, where two parameter 
estimates are obtained every time. The sample means 
and sigmas of the ESPRIT estimates were 23.93° + 
1.07° and 28.06° ± 1.37°, while those of the 3175 
successful MUSIC triáis were 24.35° + 0.28° and 
27.48° + 0.38°, Note that with reference to Fig 4, 
there is an overlap in the distributions of the ESPRIT 
estimates. This has an effect on the statistics 
calculated, since a simple angle-ordering scheme was 
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ilseü wñerein the larger of the two angle estimates in 
each trial was associated with the 28° source. The 
effect is presumed to be small in this case. 

The results indícate the presence of a bias 
even in the successful conventional MUSIC 
estimates, the source of which is described in detail 
in [ 1 ]. On the other hand, the ESPRIT estimates are 
unbiased, although of larger variance since less 
information concerning the array geometry is being 
utilized. Note also that in comparing the estimate 
variances, there is no attempt to account for the 1825 
triáis in which MUSIC failed 25 to provide two DOA 
estimates! However, as the subarray separation 
increases, the ESPRIT parameter estimate variances 
approach those of MUSIC. The same experiment was 
run for a subarray separation of 4X, and the resulting 
ESPRIT estimates were 24.003° + 0.062° and 28.002° 
+ 0.089°, and the corresponding MUSIC estimates 
were 24.011° + 0.056° and 27.986° ± 0.078°. Due to 
the increased subarray spacing (array aperture), there 
were no MUSIC failures. Again. ESPRIT is unbiased, 
but now the sample parameter estimate sigmas are 
nearly equal to those obtained with MUSIC. More 
details concerning these and many other simulations 



be 



found 



[1] 



(see 



also 
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[18]). 

VI Discussion 

In this paper, a new technique (ESPRIT) for 
signal parameter estimation has been introduced. The 
algorithm differs from its predecessors in that a total 
least-squares rather than a standard least-squares (LS) 
criterion is employed. The earlier versions of ESPRIT 
described in [15] and [16] can be seen [1] to be least- 
squares estimators of an m x m operator whose action 
is restricted to a d-dimensional subspace. The fact 
that this LS operator is a restricted m x m operator 
leads to some concern over potential numerical 
difficulties in solving the generalized eigenproblem. 
Imposing the subspace restriction as a constraint prior 
to solving the generalized eigenproblem leads to a 
well-conditioned d x d generalized eigenproblem, 
thereby mitigating these numerical concerns; 
however, the least-squares property of the estimate is 
retained. In cases where the SNR is sufficiently large, 
the difference between the LS and TLS parameter 
estimates is small. The difference is notable at low 
SNRs, however; the LS estimates are biased as 
predicted, while the TLS estimates are relatively 
unbiased and have been shown 18] to be strongly 
consistent (converge with probability one to the true 
valúes). 
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A. ESPRIT and Null-Steering 

Many of the previous high-resolution 
parameter estimation techniques are based on steering 
beams toward signal directions and. in some cases, 
simultaneously attempting to otherwise minimize the 
power in some weighted combination of sensor 
outputs. Parameter estimates are associated with 
DOA's at which peaks in the power output occur. 
Although intuitively appealing at first glance, there is 
a mueh more powerful alternative philosophy, that of 
null-steering. It is well known that deep, sharp 
notches in directivity patterns and filter gain func- 
tions are much easier to achieve than sharp peaks. 
Interferometers exploit this fact to obtain accurate 
estimates of source parameters by finding the relative 
phase required to cancel signal components in two 
channels. In this context. ESPRIT can be interpreted 
as a multidimensional null-steering parameter 
estimation algorithm. Calculation of the eigenvalues 
of the (rotation) operator i|/ which are the roots of its 
characteristic polynomial, can be interpreted as 
multidimensional null-steenng. Instead of steering 
broad beams, ESPRIT steers sharp nulls at all sources 
simultaneously and does so without relying on 
knowledge of the array manifold! 

B Computational Advantages of ESPRIT 

The primary computational advantage of 
ESPRIT is that it eliminates the search procedure 
inherent in all previous methods (ML, ME, MUSIC). 
ESPRIT produces signal parameter estimates directly 
in terms of (generalized) eigenvalues. As noted 
previously, this involves computations of the order 
d 3 . On the other hand, MUSIC and the other high- 
resolution techniques require a search over and it is 
this search that is computationally expensive. The 
significant computational advantage of ESPRIT 
becomes even more pronounced in multidimensional 
parameter estimation where the computational load 
grows linearly with dimensión in ESPRIT, while that 
of MUSIC grows exponentially. If r e is the resolution 
(i.e., number of vectors) required in the calibration of 
a for the L/th dimensión in 9, the computation 
required to search over L dimensions for d parameter 
vectors is proportional to n¡ L = 1\ r,. For r ( = r. the 
computational load is r L . 

Appendix 

The Máximum Likelihood 

Estimator For the class of problems considered 
herein, the máximum likelihood estimator is simple to 
derive analytically although, in most practical real- 
time applications, computationally prohibitive. For 
nonrandom signáis in Gaussian noise with covariance 
Z„ z(t) = A(9)s(t) + n(t), the likelihood function is 
easily written [51, (1); 

Log [z(t)] = -ln [P{z(t)lz(t) = A(9)s(t) + n}] 
(28) 

*- (z(t)-A(9) s(t)f 

-i 



Z n ' [z(t) -A(9) s(t)] 



(29) 
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"The maximization of Log is over {s (t); te [ 
o, N] } x {9 G 9} and is therefore a non linear 
optimization problem. It belongs, however, to the 
class of separable non linear optimization problem. 
Golub and Pereyra [19] prove that the optimization 
can be carried out in two setps. A solution for the 
optimal s(t) is sought as a function of 9, then the 
maximization over 9 is performed. Employing this 
procedure gives "s (t) = w* (9) z (t). Where w (9) = 

X n 1 A(9)A*(e)^ n 1 A(e)[ 1 . Substituting the 

expression for s (t) back into (29) and using standard 
properties of the trace operator, 

Log(9)<x-Tr{P A . (e) R Z z ^ 30 > 

Where Pa-(B) is the oblique projection operator into 
the complement of the space spanned by A(9) (in the 
metric Z n ). Maximization of this criterion is 
equivalent to finding. 
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max 
9 



T r {P A( e)Rzz ! (31) 



as can be easily verified. Although easy to describe 
analytically, the computational burden of actually 
carrying out the multidimensional projection and 
maximization over 9 is generally prohibitive, 
resulting in the need for reasonable approximate 
solutions such as MUSIC and ESPRIT. 
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